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Epitaxial growth via chemical vapor deposition is considered to be the most promising way towards 
synthesizing large area graphene with high quality. However, it remains a big theoretical challenge 
to reveal growth kinetics with atomically energetic and large-scale spatial information included. 
Here, we propose a minimal kinetic Monte Carlo model to address such an issue on an active 
catalyst surface with graphene/substrate lattice mismatch, which facilitates us to perform large scale 
simulations of the growth kinetics over two dimensional surface with growth fronts of complex shapes. 
A geometry-determined large-scale growth mechanism is revealed, where the rate-dominating event 
is found to be Ci-attachment for concave growth-front segments and Cs-attachment for others. This 
growth mechanism leads to an interesting time-resolved growth behavior which is well consistent 
with that observed in a recent scanning tunneling microscopy experiment. 


I. INTRODUCTION 


Graphene, a material with outstanding electronic, me¬ 
chanical, thermal, and optical performance, has gained 
explosive gr owth of research int erests si nce the works pi- 
oneered by Novoselov et al.[ T] and by [Berger et al.[ 2 ] 
Among several ways to produce graphene [HIS], epitaxial 
growth via chemical vapor deposition (CVD) on metal 
surfaces is considered to be the most promising one to¬ 
wards synthesizing high-quality and large area graphene 
at relatively low cost[4]-[8]. In order to ensure controllable 
growth of high quality graphene, great efforts have been 
contributed to understand the underlying mechanism. 0 - 
[18]. Many experimental tools are able to provide the 
structure or growth information for graphene, such as the 
scanning tunneling microscopy, low-energy and photoe¬ 
mission microscopy, surface-sensitive electron diffraction, 
and various growth behaviors have been revealed. Typ¬ 
ically, graphene growth is considered to be a precipitate 
of carbon atoms dissolved in the substrate on transition- 
metals such as Ni, Co and Fe, and carbon solubility is 
a key factor for graphene growth[T5]. On Cu surface, 
due to its low solubility and relative low catalyst activ¬ 
ity, growth process is mainly determined by diffusion of 
surface carbon species [16]. For active catalyst surfaces 
such as Ir, Rh and Ru, interaction between carbon and 
substrate atoms is strong and the graphene/substrate lat¬ 
tices are not well matched, so that the grown graphene 
spreads on the surface just like a carpet [5] and forms 
moire patterns which are usually super unit cells consist¬ 
ing of hundreds of carbon atoms [131 O [lOl- Growth 
of graphene on these surfaces needs highly supersat¬ 
urated carbon monomers and the growth rate shows 
strongly nonlinear dependence on the concentration of 
monomer m- 


In spite of plenty of experimental achievements, there 
are yet rare theoretical studies which are surely de¬ 
manded to provide deep understanding of the growth 
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kinetics during the CVD process. Generally, different 
theoretical methods can be applied to study the growth- 
related behaviors at different scales. At the atomic level, 
first principle calculations can provide detailed energetic 
information of the growth events without any empiri¬ 
cal parameters. For example, monomer (Ci) attach¬ 
ment to armchair edges of graphene is energetically fa¬ 
vorable, producing grown islands with atomically sharp 
zigzag edges [ 21 ]. Density functional theory calculations 
showed that lattice mismatch between graphene and Ir 
substrate results that Ci attachment has to overcome 
a higher energetic barrier for some growth sites than 
others [ 22 ]. Although first principle calculations are able 
to provide detailed information at the atomic scales, they 
are usually too computationally expensive to study the 
growth kinetics. On the other hand, a rate theory and 
its refinement at the macroscopic level have been devel¬ 
oped to produce a quantitative account of the measured 
time-dependent carbon adatom density by assuming that 
graphene islands grow homogeneously via the attachment 
of five-atom carbon clusters (C5)[23| [24]. With opti¬ 
mized kinetic parameters, it predicts that the smallest 
stable precursor to graphene growth is an immobile is¬ 
land composed of six C 5 clusters. Note that this model 
is macroscopic and the kinetic parameters used are ex¬ 
periential, such that revealing underlying mechanisms at 
atomic level is out of its range. 

To bridge the gap between microscopic growth events 
and macroscopic growth kinetics, very recently, we have 
proposed a multiscale “standing-on-the-front” kinetic 
Monte Carlo (SOF-kMC) approach combining with first 
principle calculations to study graphene growth on Ir 
surface [22l [25]. Therein, a quasi-lD model was estab¬ 
lished by focusing on detailed growth events of carbon 
species on the growth front to understand the special 
growth behavior where the average growth rate is a 
highly nonlinear function of Ci concentration [17]. Nev¬ 
ertheless, for general growth processes in experiments, 
there are a plenty of new features beyond the quasi-ID 
growth kinetics, which are out of the scope of SOF-kMC. 
For example, when several adjacent fragments merge into 
one whole piece of grown graphene, it has to grow over 
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a vacancy island (VI) which is an island shape of empty 
region inside the flake of graphene. In addition, inter¬ 
esting time-resolved growth behaviors involving VI over 
the large scale of moire patterns have been reported 
recently [18]. Therefore, an applicable approach for full 
2D kinetics of large-scale graphene growth is very much 
desired to obtain a deep understanding of the underlying 
mechanism, which, however, still remains a big theoreti¬ 
cal challenge. 

In this paper, we establish a minimal kMC model to 
address such an issue on an active catalyst surface with 
lattice mismatch. As shown in the inset of Figj^a), the 
model is based on a two-type-site kMC lattice where one 
type atop the underlying substrate atoms (D-site) is diffi¬ 
cult for Cl attachment and the other (E-site) is easy [ 22 ]. 
and consists of only several essential growth events dom¬ 
inating the growth process while the contribution of oth¬ 
ers are implicitly considered in effective kinetic param¬ 
eters of these dominating events. With kinetic param¬ 
eters extracted conveniently from detailed atomic cal¬ 
culations, graphene growth over VI is investigated as 
an application of the minimal model. Remarkably, we 
show that the complicated growth behavior observed in 
experiments can be well reproduced. Detailed analysis 
reveals that the interesting time dependence of the VI 
area is due to the fact that the rate-dominating events 
for growth is strongly dependent on the local shape of 
the growth front: The growth event over D-sites can be 
Cl-attachment for concave segments while has to be C 5 - 
attachment for others. Such observation suggests that 
geometry of the growth front plays important and subtle 
roles during the large-scale growth process of graphene. 


II. THE MINIMAL MODEL FOR 
LARGE-SCALE GROWTH KINETICS 

Generally, a whole process of graphene growth includes 
a nucleation stage and an epitaxial growth stage after 
nucleation, both of which are challenges for theoretical 
studies. Kinetic Monte Carlo (kMC) simulation is good 
choice for the later challenge for several reasons [23] . 
Firstly, kMC goes beyond the atomistic detail and con¬ 
siders kinetic processes in a time/space coarse-grained 
manner, which facilitates simulations of growth kinet¬ 
ics with atomic details. Secondly, wealth of spatial in¬ 
formation is contained in snapshots of kMC simulation, 
which can be directly related and compared to the ex¬ 
perimental scanning tunneling microscopy images. Nev¬ 
ertheless, as pointed out by |Tetlow et al.| building up of 
practicable kMC models for epitaxial growth of graphene 
is a remarkably complicated problem [9]. The main is¬ 
sue is about how to capture all the essential events im¬ 
portant for growth kinetics while algorithm realization 
and computational efficiency can be ensured. For large- 
scale growth of graphene, a full version of kMC model 
should include all the atomic events related to absorp¬ 
tion/desorption and dissociation of carbon sources on 
metal surface, reactions and diffusion of carbon clusters, 
attachment/detachment of carbon species to graphene 


edges, etc. Nevertheless, such a full-kMC is not applica¬ 
ble due to the disparate rates of different kMC events and 
the entanglement of detailed growth events and the large 
spatial scale of graphene. In order to overcome the first 
problem, we have built up the multiscale SOF-kMC as 
described in Ref. [22]and[25l To make the second prob¬ 
lem tractable, one should identify the very key events 
that are most relevant to the large-scale growth kinetics. 

As mentioned above, there are essentially two types of 
sites over an active catalyst surface such as Ir as shown 
in the inset of Figj^a), i.e.. E-site and D-site, as a conse¬ 
quence of lattice mismatch. As a prior step, one should 
identify possible elementary events for kMC simulations 
of large scale growth kinetics. By dividing the whole sur¬ 
face lattice into four regions: The grown graphene sheet, 
growth front, diffusion layer, and far field, the SOE-kMC 
is adopted to this end with all the attachment and de¬ 
tachment events of carbon species from Ci up to Cq while 
all other events in the diffusion layer and far field are 
compacted into effective carbon fluxes [ 22 ] l25] . An im¬ 
portant observation is that, while the growth over E-sites 
is mainly by Ci attachment, that over D-sites is a little 
more complicated |25|. As shown in the inset of Eigj^a), 
a moire pattern consists of two classes of D-sites along a 
given zigzag growth front. Eor a configuration similar to 
the one enclosed by the blue rectangle, the two E-sites ad¬ 
jacent to the completed zigzag graphene edge have been 
occupied by Ci very fast, then, attachment of Ci on the 
D-site will locally close a zigzag graphene edge which is 
stable and hard to be detached. Thus, such a D-site can 
be named as a closed D-site. On contrary, for the bottom- 
left configuration surrounded by the green rectangle, D- 
site is adjacent directly to the completed zigzag grown 
front, where the attachment of Ci will result in an open 
zigzag growth front which is very unstable and easy to be 
detached. Consequently, this D-site is denoted as an open 
D-site. Due to these observations, the overall growth rate 
is found to be dominated by and contributed slightly 
by Cq attachment on open D-sites, resulting in a highly 
nonlinear dependence of growth rate on Ci concentration 
with an exponent slightly bigger than 5 [25], while Ci at¬ 
tachment on closed D-sites also provides indispensable 
net contributions to grown graphene. 

Based on the above picture, it is convenient for us to re¬ 
fine all the attachment and detachment of carbon species 
on the growth front to be three key events: Ci attach¬ 
ment on E-sites, Ci attachment on closed D-sites, and 
C 5 attachment on open D-sites. This leads to a minimal 
model given by 


Cl ^ G 

E 

Cl > G 

^closed 


C5 



^ G 


( 1 ) 


where E, Dopen and Ddosed denote E-sites, open D-sites 
and closed D-sites, respectively, and G refers to the grown 
graphene. 
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(a) 



t (in arbitrary unit) 


FIG. 1: (Color online) (a) A typical snapshot during the sim¬ 
ulation of large-scale graphene growth with edges straight at 
the moire pattern scale. Each dot represents a moire pat¬ 
tern consisting of 10 x 20 growth sites. Black/blue dots refer 
to fully/partially occupied moire patterns, respectively. The 
inset shows a zoom in of a moire pattern, where each cir¬ 
cle refers to a growth site. The substrate atoms are in light 
gray, and E-sites and D-sites are colored by black and red, 
while solid and open cycles denote occupied and empty sites. 
The blue/green rectangle encloses a growth fragment with 
closed/open D-sites, respectively, (b) Linear dependence of 
area A (normalized by its initial value) for monomer concen¬ 
tration ni = 0.01 monolayer obtained by using the minimal 
model (Eq. 0 ). The inset presents the highly nonlinear de¬ 
pendence of average growth rate (dA/dt) on ni, which can be 
well fitted by {dA/dt) = a + bnj with an exponent 7 ~ 5. 


To accomplish our model, we need to calculate 
the effective rate constants and kc. Each 

ki{i = a^b^c) contains three components, i.e. ki ^ 
niexp{—ei/{kBT)}fo^i^ wherein rii is the surface concen¬ 
tration of the species-i, is the associated energetic 
parameters with contributions from the diffusion and 
attachment processes. Notice that the minimal model 
only considers the three dominant attaching events as 
shown in Eq. 0 , in real growth process nevertheless, 
other events should also take effects on the growth pro¬ 
cess. For example, the attached carbon atoms may be de¬ 
tached via detaching events, which leads to a decreasing 
of the effective attaching rate. To take this into account, 
a prefactor /o,i denoting a net contribution of event i 
to the front growth is added as the third component of 


the effective growth rate constant. As already described 
in detail in our previous studies [ 22 l |25], the energetic 
parameter associated with can be calculated by first 
principle calculations and rii can be derived by assum¬ 
ing a quasi-equilibrium between differently sized carbon 
species. The factor /o,i can be extracted from detailed 
simulations with full attachment and detachment events 
on a small sized lattice as shown in Refl25l 

Since all other events, such as combining, diffusion and 
detachment of carbon species, have been compacted in 
the effective rate constant ki^ a standard kMC simula¬ 
tion is sufficient for the simulation of large scale growth 
of graphene. One only needs to figure out configurations 
on the growth front ready for attachment, and the rate 
for each event is = kiNi where Ni is the number of 
configuration for event i. To check the validity of our 
approach, we first apply the above minimal model to in¬ 
vestigate graphene growth with steady carbon fluxes and 
growth front that is straight at the moire pattern scale on 
Ir(lll) surface. As illustrated in Figj^a), kMC is per¬ 
formed on a huge lattice containing x = 80 x 48 
moire patterns each of which consists of 10 x 20 growth 
sites with an initial grown graphene ribbon of size 2 * 48 
moire patterns. The rate constants are ka = 2.18 x 10^, 
kb = 2.18 and kc = 2.455 x 10“^ for monomer carbon 
concentration rii = 0.01 monolayer and T = llTOiF by 
taking /o,c to be the unit which leads to the magnitude 
of /o,a is about 10^ and /o ,6 ~ 1. In practice, we find 
that the growth kinetics is very robust to the exact val¬ 
ues of these rate constants due to the large discrepancies 
between them. 

In this case, the rate-dominating step is the attach¬ 
ment of C 5 on the D-sites of open-type. Therefore, the 
growth rate of the graphene area dA/dt is proportional 
to the number of configurations on the growth front that 
allow the rate-dominating event happening (denoted here 
by Urdc)- In Fig 0 b), we plot A (normalized by the 
area of initial grown ribbon) as a function of time t, 
where a perfect linear dependence is observed, indicating 
Urdc are nearly constant along the growth front. This 
is reasonable since Urdc should be proportional to the 
length of the growth front, here is just Ly which is a 
constant. In the inset of Figj^b), time averaged area- 
growth rate {dA/dt) is plot as a function of the monomer 
concentration ni. The curve obtained can be fitted by 
{dA/dt) = a bn/ with an exponent 7 ~ 5, which is 
well consistent with observations in experiments pT] and 
in previous simulation [2^. Thus we believe that the 
minimal model is able to capture the essential physical 
picture regarding the growth kinetics on Ir surface with 
graphene-substrate lattice mismatch. 


III. APPLICATION ON GROWTH OVER 
VACANCY ISLAND 

Now we use our approach to investigate the growth 
behavior with complex front shapes. As a typical exam¬ 
ple, we consider a VI as depicted in the inset of Figj^a) 
wherein several fragments merge into on piece of grown 
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which gives 

+ nt (3) 

where, m, n and p are certain fitting parameters. Re¬ 
markably, we find that such a function can fit the ex¬ 
perimental data quantitatively well as demonstrated in 
FigJI^a) with m = -3.74 x 10“®, n = -5.47 x IQ-^ 
and p = 1.02. Therefore, we believe that Eq.([^ and § 
have captured some basic physical features regarding the 
growth of graphene over the VI. 





FIG. 2: (Color online) (a) The area A obtained experimen¬ 
tally by |Dong and Frenken[l8j (black squares) as a function 
of time t for graphene growth over VI, which can be well fit¬ 
ted by Eq. § (the red line) derived from the minimal kMC 
simulation. The inset shows the initial shape of VI at the 
moire pattern scale, (b) Linear dependence of the growth 
rate dAvi/dt normalized by y/Avi on t. The inset shows the 
simulated time-resolved Ayi of VI. 

graphene. The growth kinetics over such a specific shape 
has been investigated in some details experimentally [18]. 
It has been found that the area Ayi (the area of the 
empty region inside the graphene flake) normalized by 
its initial value exhibits a complicated dependence on 
time, which is sur ely non-exponential as discussed by 
Dong and Frenk^lSj . while the underlying mechanism 
is yet not clear. Interestingly, using our approach with 
all the parameters same as above mentioned, we can ob¬ 
tained the time-dependence of Ayi (the inset of FigJ^b)) 
which is very similar in tendency with that observed 
in the experiments (scatters in Figj^a)). In addition, 
our model makes it convenient to investigate in detail 
how the boundary of the VI change with time. As dis¬ 
cussed above, for a regular front as shown in Figj^a), one 
can expect that dAyj/dt oc Ayj since rirdc oc y^Ayj. 
Nevertheless, for irregular boundary here, we And that 
dAyj/dt is not proportional to \/Ayi at all. Rather sur¬ 
prisingly, {dAyi Jdt) / \/ Ay I in nearly linear with time t, 
as drawn in FigMb). One may then write 

= yj Avi{t){mt + n), (2) 


FIG. 3: (Golor online) Typical kMG snapshots for growth over 
VI at (a) t = 17.8 (in arbitrary unit) and (b) t = 40.7 at the 
moire pattern scale. Red lines illustrate the initial shape of 
VI. (c) and (d) are typical zoom-in atomic-level structures for 
a straight segment and a concave one, respectively, where each 
circle refers to a growth site and E-sites/D-sites is colored by 
black/red. The green rectangles enclose typical configurations 
of rate-dominating events for growth over D-sites. In (d), 
D-site is actually of closed-type along a given zigzag front 
indicated by Ri, although along R 2 it is of open-type. 

To obtain a deeper understanding about the interesting 
growth behavior of VI, we try to And hints from spatial 
information contained in snapshots of the surface dur¬ 
ing growth process, which is one of the major advantages 
of kMC simulations. Two snapshots of VI are shown 
in FigJ^a) and (b) for t = 17.8 and t = 40.7 (in arbi¬ 
trary unit), respectively, where the initial shape is out¬ 
lined by the red lines. Compared to the initial profile, 
it can be observed that the growth process mainly takes 
place at the concave segments in the corner parts, while 
the convex parts in the middle remain nearly unchanged. 
This indicates that the number of relevant configuration 
Urdc should be proportional to the length of concave seg¬ 
ments within the growth front, Lc? rather than to the 
total length of growth front L. Comparing the snapshots 
for t = 17.8 and t = 40.7, one can see that the ratio L^jL 
increases with time. To be in accordance with EqJ^, one 
expects that L^jL oc t assuming that L oc Ay ^ Nev¬ 
ertheless, the reason of this particular linear dependence 
with time is still open to us at the current stage. 

Now the key point turns to why the effective num¬ 
ber of configurations Urdc for growth is proportional to 
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the length of concave front segments. To this end, we de¬ 
pict typical zoom-in atomic-level structures for a straight 
segment and a concave one in FigJ^c) and (d), respec¬ 
tively. Clearly, for the straight segment, the D-site is of 
open-type, such that the rate-dominating step still has 
to be C 5 -attachment. However, for the concave segment 
formed within the corner environment, the D-site is actu¬ 
ally of closed-type along a given zigzag front indicated by 
i?i, although along R 2 it is of open-type. Such a dynamic 
change from open-type to closed-type of the D-site due to 
the local shape of the front leads to a change of the rate- 
dominating event from Cs-attachment to Ci-attachment. 
In other words, the rate-dominating event (growth over 
D-sites) is heterogeneous for different parts of the growth 
front, i.e., Ci attachment for concave segment, and C 5 
attachment for others. Consequently, the overall growth 
rate of graphene is determined by the length of concave 
segments. Noted that the energetic parameters for all 
growth events are independent on the front shape, the 
heterogeneity of rate-dominating events is solely geome¬ 
try effects. 

IV. CONCLUSION 

In summary, we have built up a minimal kinetic 
Monte Carlo model to study large-scale growth kinet¬ 


ics of graphene on active catalyst metal surfaces with 
graphene/substrate lattice mismatch. Using kinetic pa¬ 
rameters extracted from detailed atomic calculations, 
our model successfully reproduced the time-resolved 
graphene growth over vacancy island observed in ex¬ 
periments, and revealed that the interesting dependence 
of the area on time is resulted from the heterogeneous 
rate-dominating events for graphene growth at different 
growth front segments, which was shown to be mainly 
determined by geometry. Since the proposed minimal 
kMC model (or perhaps its modified version) provides a 
powerful theoretical tool for investigation of large-scale 
growth kinetics of graphene, and our findings take an 
important step forward to the underlying growth mecha¬ 
nism, we hope that the present study could open new 
perspectives on theoretical studies and motivate more 
atomic-resolution imaging experiments. 
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